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Abstract 

Continuing our investigation into the Hierarchical Reference Theory of 
fluids for thermodynamic states of infinite isothermal compressibility kt 
we now turn to the available numerical evidence to elucidate the character 
of the partial differential equation: Of the three scenarios identified previ- 
ously, only the assumption of the equations turning stiff when building up 
the divergence of kt allows for a satisfactory interpretation of the data. 
In addition to the asymptotic regime where the arguments of part I di- 
rectly apply, a similar mechanism is identified that gives rise to transient 
stiffness at intermediate cutoff for low enough temperature. Heuristic ar- 
guments point to a connection between the form of the Fourier transform 
of the perturbational part of the interaction potential and the cutoff where 
finite difference approximations of the differential equation cease to be ap- 
plicable, and they highlight the rather special standing of the hard-core 
Yukawa potential as regards the severity of the computational difficulties. 

Keywords: liquid- vapor transitions, non-linear partial differential equations, nu- 
merical analysis, finite differences, stiffness. 

1 Introduction 

Building upon the results of part I [1], q. v., and maintaining the notational and 
semantic conventions introduced there, we now turn to the numerical solution 
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Table 1: Overview of systems and sample isotherms: (3 C and g c give the location 
of the critical point, g v and gi the extent of the two-phase region at the inverse 
temperature [3 considered in the tables and figures to follow. The numbers have 
been obtained from hrt calculations not imposing the core condition. All of 
the digits indicated for f3 c are significant. 
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of the hrt pde [2, 3, 4, 5, 6, 7, 

df , , 8 2 f 
&Q- doo + do2 di 2 

in order to determine the type of behavior that actually occurs in practical 
applications of the theory for thermodynamic states of diverging isothermal 
compressibility kt- To this end we consider two simple model potentials v(r) = 
v hs (r) +w(r), viz., the hard-core Yukawa (hcy) system, 



and square wells (sws) 



■ e -^'-°> : r > a, 



to illustrate the types of behavior encountered and to test the predictions fur- 
nished by the relevant scenarios. In both of these potentials e coincides with 
the negative of the contact value of the interaction, lim T _ y(T+ (— w(r)), and so 
sets the energy scale of the problem. The potential range is given by 1/z and 
Act, respectively. Unless stated otherwise, eo, the value of w hcy (r) inside the 
core, coincides with e, a choice shared with the implementation by the authors 
of HRT and their coworkers referred to as the original one in refs. [9, 10], q. v.. 
A short summary of the parameter sets and sample isotherms considered in 
this study can be found in tab. 1. In the numerical work we employ an un- 
conditionally stable implicit predictor-corrector scheme shortly characterized in 
section 3.1. A more extensive discussion of the implementation can be found 
in refs. [9, 10], where default settings for the most important customization pa- 
rameters arc also documented. Even further technical information is available 
with the source distribution itself [11]. 

Of the three types of behavior compatible with the local properties of the 
PDE, both genuine (r = s = 0) and effective (r > 0, s = r+1 > 1, r off = s off = 0) 
smoothness imply a fd approximation to / growing like 1/Q. The monotonous 
scenario, on the other hand, furnishes the specific prediction that E Q 2 tends to 
a finite limit for Q — > 0. As we will see in section 2, the numerical evidence 
clearly rules out this possibility. 
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It is thus only the genuinely smooth and the stiff scenarios that remain to 
be considered in section 3. The results of the computations reported there do, 
indeed, allow us to infer the character of the pde for subcritical temperatures, 
T < T c , with great confidence, if only indirectly due to the great computational 
similarity of genuine and effective smoothness. Our main evidence in favor 
of the stiff scenario derives from the rather detailed and testable predictions it 
entails, all of which are confirmed numerically. By way of contrast, the genuinely 
smooth scenario does not hold an explanation for the observed trends, especially 
as regards the dependence of the fd results on the properties of the discretization 
grids. 

Our conclusion that the PDE actually turns stiff in part of V for T <T C then 
paves the way for some heuristic arguments relating the onset of smoothing in Q 
to the form of the Fourier transform of the perturbational part of the potential 
(section 4). So having understood the behavior of the pde in the limit Q — > 
where asymptotic reasoning valid for large e applies, in section 5 we then turn 
to similar computationally problematic features of its solution at much higher 
cutoff where the numerical evidence points to a mechanism not unlike that at 
work in the asymptotic region. We close with an informal discussion of the 
reasons for the atypical computational properties of HCY fluids of moderate 
inverse screening length z (section 6). 

2 The monotonicity assumption refuted 

According to section V of part I [1], the assumption of a merely logarithmic 
divergence of / furnishes the rather specific prediction of eQ 2 tending to a 
finite limit for Q — > 0. Of course, the possiblity of non-zero s means that, in 
principle, the smoothing effect discussed in section VI of part I [I] must be 
reckoned with. The singularity being so mild, however, a possible reduction of 
s > to an effective value of s e s = is preempted by the choice of step sizes 
AQ: 

In our implementation of the theory the cutoff in the i th FD step is parametrized 

as 

Q (l) =\n(e a - lb + 1) /a, i = 0, 1, . . . , 

as is the case for the program the original authors of hrt and their coworkers 
employ, too. Here aja is close to the cutoff Qr \ = Q x where initial conditions 
are imposed on /, and b/a is the spacing of successive cutoffs in the large 

Q limit. For Q — > 0, on the other hand, we easily find AQ ~ —Q (f — e~ b ). If 
e Q 2 is to approach zero or a finite constant as predicted by the assumption of 
monotonous growth these step sizes thus turn out of order 0(e -1 / 2 ) at most, 
and our discretization should allow us to follow the variation of / reasonably 
well all the way to Q = 0. From fig. I, however, we see that eQ 2 clearly 
diverges for Q — > 0. As this finding is corroborated by further calculations 
with a smaller setting of the numerical parameter AQ | , on finer density grids 
(down to Ag = 5 • 10~ 4 /cr 3 ), and for both hard-core Yukawa and square well 



3 




0.2 




P = 
P = 
P = 
P = 



0.2/d 3 
0.3/d 3 
OA/a 3 
0.5/d 3 



• A 

Y 



0.3 



Figure 1: eQ 2 as a function of Q for various densities inside the binodal: The 
data have been obtained for a hard-core Yukawa potential with inverse screening 
length z = 1.8/<7 and for an inverse temperature of (3 — 0.875/e. The numerical 
precision in the calculations was e# = 1CP 2 , the step size for infinite cutoff was 



potentials we feel we can safely exclude the monotonous growth scenario from 
further consideration. 

3 Smoothness vs. stiffness 

As for the remaining two alternatives, an attempt to distinguish numerically 
between genuine and effective smoothness seems doomed at first sight: both 
predict a numerically smooth solution growing like 1/Q and with a profile like 
that sketched in fig. 1 of part I [1]. And indeed, fig. 2 shows the small Q behavior 
of / within the binodal as obtained numerically to be in excellent agreement 
with / oc 1/Q, and figs. 2, 3, and 5 as well as the numerical data demonstrate 
that / is of the form necessary for a stable pattern of growth, as postulated in 
part I [1], q. v.. But even though these general features fit both scenarios, close 
scrutiny of the computational process and the numerical results yields a wealth 
of indirect evidence that we feel is sufficient to establish the stiffness of the PDE 
for T < T c with great confidence, if not with absolute certainty. We base our 
reasoning on the rather specific, and numerically testable predictions that follow 
from the assumption of stiffness. These stand in marked contrast to the vague 
expectations furnished by the genuinely smooth scenario that can, however, 
never be ruled out completely as smoothness is its sole defining characteristic. 
As we will show in this section, it is the assumption of a stiff pde that is in full 
accordance with the numerical findings whereas smoothness is only marginally 
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Figure 2: / as a function of Q for various densities inside the binodal: The 
data have been obtained for the same hard-core Yukawa potential and with 
the same numerical parameters as in fig. 1. The dashed line indicates the slope 
corresponding to proportionality of / to the reciprocal of the cutoff. Subsequent 
symbols are separated by ten steps in the — Q direction. 



compatible with some of their traits, especially as regards the SW data. 

Before, however, some general remarks are in place: Letting the labels x and 
y refer to either Q or g, in the stiff scenario smoothing in x sets in at Q = Q^x 
and can always be postponed, i. e., shifted to lower cutoffs by decreasing the step 
size Ax. If, however, the corresponding exponent, r or s, is positive, the rapid 
growth of /, e and \d 2 f /dx 2 \ as the solution proceeds towards Q = implies 
that the amount by which Qax can be changed in this way and the attendant 
computational effects must be small. For positive exponents the Qax are thus 
fairly well defined despite the gradual nature of the transition to the smoothing 
regime. — Furthermore, without loss of generality assuming Qax > Qa v , Qax 
is obviously independent of the step sizes Ay. The solution obtained numerically 
at cutoffs below Q sm0 oth = Qax is already affected by smoothing in x so that 
there is no point in identifying Qa v with the cutoff where Ay becomes too large 
to describe the variation of the no longer accessible true solution of the pde. 
Instead, Q\ y is taken to be the cutoff where smoothing in y commences in the 
solution of the fd equations (fdes) , which implies a Ax dependence of Qa v and 
may even induce Qa v to vanish altogether. For Q < Qa v , the solution generated 
numerically by necessity conforms to the smooth scenario as r e s — s cff = and 
so grows like 1/Q in a stable manner. This proportionality also means that the 
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form of / remains constant from QAy all the way to the final results at zero 
cutoff. (Here and below the form of / at some cutoff Q refers to f(Q, g) as 
a function of g, restricted to gi < g < g 2 and without regard for the overall 
normalization of /.) — The important mechanism sketched in section III of 
part I [1] and the concomitant stabilization of form and monotonicity of / do 
not explicitly depend on s and thus always set in at Qa b ', incidentally figs. 2 
and 3 show its preconditions, viz., flatness and compatibility with the sketch 
of part I [1] to be met numerically. Of course, both Qaq and Qaq depend on 
temperature and density, which is taken to be silently accounted for whenever 
we speak of the form of / at one of the Qax, and they are defined only in that 
part of V where / is large. — Not surprisingly, the two possible orderings for the 
cutoffs Qaq and Qa b assigned in an interpretation of the numerical results in 
terms of the stiff scenario entail vastly different consequences and are therefore 
discussed separately in subsections 3.2 (Qa b > Qaq) and 3.3 (Qaq > Qa s ) 
below. 

Before that, however, it is worthwhile to step back for a moment and ask 
why we have to adopt eq. (1) in the first place if the most direct formulation 
of the theory is that of a PDE for the free energy A^\g) of the Q system at 
density g, cf. part I [1]. Indeed, from eqs. (A2) and (A3) of part I [1] we see that 
jdQ oc Q 2 (/ + const) for Q a <C 1 so that the Q and g scales characteristic 
of A(Q\g) are essentially the same as for f(Q, g). In the smooth scenario there 
is then no reason for the formulation in terms of / to be preferable to that 
in terms of the free energy, provided proper care is taken to ensure stability 
and convergence. This has certainly been the case in our earlier work shortly 
summarized in appendix B.l of rcf. [10] that nevertheless was unable to proceed 
to small Q for T < T c . Similar difficulties are reported in ref. [6], and to the 
best of our knowledge there are no hrt results on simple one-component fluids 
for T < T c except in the quasilinear formulation of eq. (1) or variants thereof. 
- In the stiff scenario all this is, of course, to be expected as the rapid low 
amplitude oscillations of the solution in this case necessitate step sizes that are 
reduced as some inverse power of e or the exponential of dA^ /dQ, and only 
under special circumstances do the discretized equations allow one to obtain a 
solution with the much larger step sizes used in practical applications. As noted 
in section II of part I [1], the auxiliary quantity f(Q,g) was introduced exactly 
for this reason [8]. 

3.1 Numerical aspects 

As some of the numerical effects are rather subtle, we should also recall several 
key aspects of the implementation we rely on. This is a highly flexible and 
fully modular computational framework for the solution of a FD approximation 
of the PDE by an implicit predictor-corrector scheme thoroughly discussed in 
refs. [9, 10]. For consistency with part I [1], in the calculations reported here we 
refrain from implementing the core condition. The discretization is applied on 
uniform density grids and with the predetermined step sizes AQ of section 2. 
Convergence of the FD equations has been checked, and iteration of the corrector 
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step does not bring about noticeable changes. 

In practical applications, the discretized equations generally cannot be solved 
down to arbitrarily small Q for T < T c , and the smallest cutoff reached we 
denote Q m in- As the failure modes responsible for an end of the program are 
known [9, 10] and can be linked to the local behavior of the solution, v. L, the 
systematic changes in Q m i n upon variation of aspects of the numerical procedure 
provide a powerful and readily accessible diagnostic tool. For the calculations 
reported here, the immediate cause for abortion of the computation at some 
cutoff Qmin is either an insufficient adaptation of the rescaling necessary for 
representing quantities affected by exponentiation of / — the scale of fig. 1 
alone shows that, e. g., e cannot be represented in double precision — or else 
because of non-real / and negative e = e + 1 in the predictor step. These two 
effects are linked to rapid increase and decrease of /, respectively. 

Unlike the Qax, Qmin obviously does not depend on the density. Instead, it 
is essentially determined by the physical potential w(r), the temperature, the 
discretization grid, and the formulation of the theory [9]. As for the latter, if 
the PDE is coupled to further constraints, and the solution vector augmented 
by additional components to be determined accordingly, the likelihood of an 
early termination of the computation in the predictor step generally increases, 
and so does Q m in- As the customary manner of implementing the core condition 
involves an expansion of the direct correlation function inside the core [6, 9], the 
sensitivity of Q m in to an increase in A cc , the number of expansion coefficients, 
again proves of interest. 

For a more detailed account we refer the reader to refs. [9, 10] as well as the 
documentation that comes with the source code distribution itself [11]. 

3.2 Smoothing in g first 

So let us first turn to the hcy fluid of inverse screening length z — 1.8/ a 
already considered in ref. [9] . As mentioned before, the numerical solution must 
be smooth at any rate and is therefore always compatible with the genuinely 
smooth scenario. In this case we expect only a small dependence of the results 
on AQ and Ag that should be essentially stochastic in nature, stemming from 
the truncation error in an otherwise unproblematic fd approximation of the 
pde alone. 

As we shall see in a moment, the numerics can also be reconciled with the 
stiff scenario if only we assume smoothing to occur in the g direction first, 
Qaq > Qaq, furnishing the following predictions: The mechanism responsible 
for stable growth of / (cf. section III of part I [1]) being at work at all cutoffs 
below Qsmooth, the stability of the computational process is not an issue and 
incorporation of the core condition is entirely unproblematic. An overflow due 
to an insufficient adaptation of the re-scaling of non-O(l) quantities is the only 
possibility for numerical failure, and its likelihood is greatly reduced when AQ 
is decreased so that smaller step sizes are generally accompanied by smaller 
values of Q m in- A systematic Ag dependence of Q m ; n is not anticipated. — For 
fixed density grid, Q\ e = Q smoo th cannot depend on AQ, nor can / at Qaq- 
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Table 2: AQ dependence of the final results for a hard-core Yukawa system: 
Just as in figs. 1 and 2, z = 1.8/a, = 0.875/e, and Ag = lir 2 /^. We use 
the notation f x for f(Qmin,x/a 3 ) and define F v x = f y /f x - 
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Table 3: AQ dependence of the final results for a hard-core Yukawa system: 
The parameters and notation coincide with those of tab. 2, except for Ag — 
5-10- 4 /(7 3 . 



On the other hand, even though smaller step sizes, i. e., smaller AQ| , cf. 
section 2, correspond to smaller Qaq, s > 1 implies that the drop in Qaq must 
be exceedingly small. As furthermore the evolution from Qa b down to Qaq is 
determined by the solution at the onset of smoothing and the properties of only 
the density grid, the form of / below Qaq, including Q = 0, is virtually AQ 
independent. — As for a variation of the density grid at fixed AQ, a reduction 
of Ag clearly entails a shift of Q sm0 oth = Qaq to smaller cutoffs, which may in 
turn cause a change in Qaq, too. These effects must be rather small because 
of the non-zero cxponenets r and s, and they must vary with the density for 
the same reason the Qax are density dependent. A change of the g grid thus 
implies a small change of the form of / at Qaq and, hence, at Qaq and all 
smaller cutoffs. As long as Qa b does not fall below Qaq, however, the ratio of 
the forms of / as obtained on different density grids cannot depend on AQ. 
All these predictions are confirmed in the actual calculations for a hcy po- 
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Table 4: Ag dependence of the form of the final results for a hard-core Yukawa 
system at varying AQ: The parameters coincide with those of tabs. 2 and 3. 
Perusing the notation introduced there, G x is Q m i„ a as evaluated for Ag = 
5 • 10~ 4 /cr 3 divided by the same quantity for Ag = 10~ 2 /a 3 . 
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tential with z = 1.8/c on density grids with A g = lO~ 2 /a 3 and Ag = 5T0~ 4 /er 3 
and varying AQ as summarized in tabs. 2 to 4 and fig. 2. For fixed density grid 
(tabs. 2 and 3, respectively), Q m i n and, hence, the final values of / markedly 
depend on AQ, the former generally decreasing and the latter increasing upon 
reduction of the step size. On the other hand, both the form of / and its mag- 
nitude at fixed cutoff — to be found in the tables under the headings of F v x and 
fx Qmin c, respectively — remain largely unchanged. Comparing the results ob- 
tained with different settings for Ag, the change in the final values of / is indeed 
almost completely due to the differences in Q m i n . The magnitude at fixed Q, on 
the other hand, is affected only moderately, viz., by a few per cent for a twenty- 
fold increase in the density resolution, and it depends on g but not on AQ, cf. 
tab. 4. Q m i n itself is not affected by the density grid in a systematic way. There 
are two sample isotherms, viz., the ones at AQ^ = 0.003/(T 3 , Ag = 10~ 2 /er 3 
and at AQj^ = 0.004/(T 3 , Ag = 5-10~ 4 /er 3 , that founder at comparatively large 
cutoff. Of these, only the former does not enter the asymptotic regime where 
/ oc 1/Q, as can clearly be seen from tab. 4. All in all, the numerical results are 
in excellent agreement with stiffness, and we note that for this system and the 
density grids considered Qaq must be sought around 10 _2 /cr. Trivially, being 
smooth the results also conform to the smooth scenario as mentioned before. 

3.3 Smoothing in Q first 

In our previous work on HRT [9, 12] we repeatedly stressed the vastly different 
numerical properties of the hcy and sw potentials. This is certainly not an- 
ticipated for genuine smoothness that merely predicts a small grid dependence 
stemming from the local truncation error of the discretization, exactly as for 
the hcy system. Still, the assumption of a genuinely smooth solution is cer- 
tainly compatible with the numerics, if only marginally so in the face of the 
most prominent feature of the evolution of /, viz., episodes of much more rapid 
variation than mere proportionality to 1/Q. 

Assuming the pde to turn stiff for large / instead, and furthermore Qaq 
to exceed Qaq for the present system, there is a Q range Qaq < Q < Qaq 
where fdes are used with inappropriately large step sizes AQ while oscillations 
in g arc not yet suppressed. For these cutoffs, the stabilization wrought by 
the mechanism analyzed in section III of part I [1] is not effective yet, and 
there is no reason for / to be convex from below throughout the density range 
gi < g < Q2- On the other hand, the overall profile of / is expected to resemble 
fig. 1 of part I [1] , and s e s = once more suggests a general growth proportional 
to 1/Q. The sign of d 2 f/dg 2 is thus unconstrained, and its modulus increases 
in unison with /, i. e., in proportion to 1/Q. As d^/d^ is of order O(l), 
however, the O(l) growth of d 2 f/dg 2 may well be sufficient to destabilize the 
growth at some Q e {Qa b ,Qaq), prompting much more rapid variation of / 
as a function of Q. Of course, these near-discontinuities of / will occur at 
different cutoffs for different densities, most often close to g\ and Q2 where the 
Qax are smallest, and neighboring densities will experience them at roughly the 
same cutoff. Furthermore, in principle the jumps should lead to both increases 
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Table 5: AQ dependence of the final results for a SW system: Just as in fig. 3, 
A = 3, [3 = 0.115/e, and Ag = 10~ 2 /er 3 . The notation is the same as in tab. 2. 
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Table 6: AQ dependence of the final results for a SW system: The parameters 
and notation coincide with those of tab. 5, except for Ag = 5 • 10~ 4 /(7 3 . 

and decreases in /, depending on the sign of d 2 f/dg 2 at slightly larger Q. 
Considering the numerics, however, a large change in / is almost certain to bring 
the calculation to an end, and all the failure modes discussed in section 3.1 are 
relevant for Q m ; n . A comparatively mild increase of /, on the other hand, may 
relax the relative curvature of / to the point of allowing the solution to enter 
once more an episode of near-stability characterized by growth in approximate 
proportion to 1/Q. As for an incorporation of the core condition, in accordance 
with section 3.1 the attendant introduction of additional degrees of freedom 
is likely to exacerbate the risk of triggering such a jump in /, c/. section 5. 

To understand the grid dependence of the numerics under the assumption 
of stiffness, recall that Q m i n itself is the location of a failed jump in /. As 
smoothing in Q is the driving force behind the computational process, Q m i n 
must be quite sensitive to AQ, but there is no reason for Q m i n to be monotonous 
in AQ. The density grid, on the other hand, is still adequate for the elliptic 
boundary value problem in g at constant Q. If the numerical process were 
stable, there should thus be no appreciable dependence of the results on Ag at 
all. In the absence of the stabilization wrought by smoothing in g, however, 
even the small differences seen upon variation of Ag must be expected to shift 
the episodes of rapid evolution to slightly different cutoffs in an unsystematic 
way. By the same token, the Ag dependence of the final form of / should be 
small, and different AQ should leave it unaltered as long as the number and the 
approximate positions of the jumps do not change. As those are least frequent 
close to the maximum of /, its form is expected to be most stable in the central 
part of the density interval of large /. 

Again, these predictions compare favorably with the numerical results for a 
SW potential of range A = 3 obtained on the same discretization grids as the 
HCY data of section 3.2. The most prominent feature, barely compabitle with 
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Figure 3: / as a function of Q for various densities inside the binodal: The 
data have been obtained for a SW potential with A = 3 and at a temperature of 
(3 = 0.115/e; otherwise, the remarks of fig. 2 apply. Arrows mark several of the 
near-discontinuities discussed in section 3.3. 

genuine smoothness, viz., near-discontinuities of / can actually be found in the 
numerical data underlying tabs. 5 and 6 at the locations marked with arrows 
in fig. 3; indeed, several of them can be seen clearly even on the logarithmic 
scale of the graph. All the other consequences of stiffness with Qaq > Qa 6 are 
also in agreement with the data of tabs. 5 and 6: In particular, a pronounced 
AQ dependence of Q m i n is accompanied by only a very modest effect as Ag is 
varied, even though the relative change in Ag is much larger than that in AQ. 
Excluding the pathological data with negative / (v. £.), the final forms of / 
are mostly AQ independent, and the forms obtained on the two density grids 
differ but slightly. Only the isotherm with AQ\ x = 0.010/er in tab. 6 presents 
a somewhat different shape than those at smaller AQ| . The differences in 
the numbers given under the heading F% are, however, still in accordance with 
the stiff scenario as discrepancies appear only close to the edge of the density 
range of large /. As for the first entry of tab. 5 (AQ| = 0.003/cr), negative 
/ corresponds to exceedingly small values of e = e + 1 ~ 10~ 27 . This is found 
to be the result of a downward jump from / ~ +10 4 (e ~ 10 5000 ) at only 
slightly higher cutoff where the form of / again corresponds to that of the other 
isotherms. Clearly, even a minor perturbation of the numerical process might 
easily have led to negative e and hence to a numerical exception; in this case our 
implementation would have discarded the last step, and the final results would 
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once more conform with those of the remainder of tab. 5. 

Let us shortly return once more to the most salient feature of the numerical 
solution, viz., its near-discontinuities. Disregarding the analytical considerations 
of part I [1] it might be tempting to imagine that, for T < T c , the pde generates 
a shock front approximately symmetrically moving outward towards the densi- 
ties q v and gi of the coexisting phases as Q approaches zero. In this view of the 
numerical process the jumps occur when the shock reaches the corresponding 
density. Such an interpretation is not consistent with the data: According to 
fig. 3 the near-discontinuities of / occur repeatedly at the same density (most 
conspicuously for g = 0.1/er 3 ), and rapid change at one density is generally ac- 
companied by similar behavior at other densities. Neither of these observations 
is compatible with the idea of a moving shock front, nor is there any reason why 
the binodal should be linked to a shock front in SWs but not in the hcy fluid, 
cf. section 3.2. 

3.4 Assertion of stiffness 

Summarizing the numerical evidence presented so far we find that of the three 
scenarios found in part I [1] only the possibility of a merely logarithmic singu- 
larity of / can be ruled out with certainty. We are then faced with the two 
alternatives of genuine smoothness of the pde on the one hand, and effective 
smoothness as a result of an fd approximation to a stiff PDE on the other hand. 
As shown in the preceding subsections 3.2 and 3.3, neither of them is in direct 
contradiction with the numerical data. 

The crucial difference is their respective specificity and testability: The gen- 
uinely smooth scenario does not make any predictions beyond the smallness of 
the discretization grid dependence of the numerical results, nor does it offer any 
of the detailed understanding of the computational process that is necessary 
for accurate and reliable interpretation of the FD results. By way of contrast, 
stiffness of the pde in part of its domain provides a consistent framework for 
the interpretation of the numerics and furthermore entails a number of concrete 
and numerically testable consequences, all of which are in excellent agreement 
with our data once the correct ordering of the Qax has been chosen. In combi- 
nation with the analytical considerations of part I [1] and our earlier statements 
regarding the importance of the formulation of the hrt pde employed, the speci- 
ficity and great number of these predictions provide ample, although necessarily 
indirect evidence in favor of the stiff scenario. 

From this point on we will therefore take it for granted that the HRT pde 
does, indeed, turn stiff in part of its domain for subcritical temperatures. On 
this basis we now aim to further enhance our understanding of the hrt numerics, 
shedding some light on the location of Qaq (section 4), extending our findings 
in the asymptotic region to intermediate Q (section 5), and finally clarifying the 
outstanding numerical properties of the hcy potential vis-a-vis other physical 
systems (section 6). 
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4 The onset of smoothing in Q 



Considering the great importance of the relative order of the Qax for the numer- 
ical process, it is natural to inquire into their typical values. As the exponents r 
and s are non-zero by assumption, these cutoffs may only weakly depend on the 
discretization grid and so are largely determined by the perturbational part of 
the potential alone. Whereas the onset of smoothing in g eludes simple reason- 
ing so far, some heuristic argucmnts point to a simple connection between the 
likelihood of finding Qaq at some cutoff and the form of the Fourier transform 
w(k): 

Let us consider a thermodynamic state of diverging isothermal compressibil- 
ity at a cutoff that is low enough for smoothing in Q to have set in at least par- 
tially, Q ~ Qaq- In view of the gradual transition between the smoothing and 
non-smoothing regimes, the effective exponent s e a may not vanish exactly yet; 
nevertheless it seems safe to assume s c ff < 1 ■ Of course we expect e » / ^> 1 so 
that reasoning based on the asymptotic behavior for large e is applicable, and 
due to the monotonicity of the exponential function the likelihood of finding 
Qaq close to some cutoff Q increases with the slope —df/dQ of /. At the same 
time, for a hard-sphere reference system Qaq can only depend on the form of 
the Fourier transform of the perturbational part of the interaction potential, 
i. e., on uq = 4>/<fio rather than on <j> itself: The temperature T = l/fcs/3 enters 
the calculation only as a pre- factor to the interaction potential, viz., through 
<f> = —f3w so that the normalization of <p only fixes an energy or temperature 
scale. 

With this in mind we define an auxiliary quantity ip(Q, g), corresponding to 
4>o + 7o^' m the notation of our earlier work on hrt [1, 9, 10, 12], through 



Solving this relation for ip and differentiating with respect to Q we obtain 



which is valid at all cutoffs except close to the zeros Q^ ti oi(j) and uq where eq. (4) 
cannot be inverted. — An alternative expression for dip/dQ can be obtained 
from the pde (1) and the compressibility sum rule: Following section 2.4.1 of 
ref. [10], for density independent potential we easily find 



Equating these two expressions for dtp/dQ, solving for d 2 f/dg 2 , and inserting 



JC + lpUQ = --. 



(4) 




dQ ~ 4?9^ 
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Figure 4: m as a function of Q for sw and hcy potentials with the parameters 
indicated: If, contrary to eq. (2), the Yukawa form is retained even inside the 
core, u (Q) is given by z 2 / (z 2 + Q 2 ). As far as the SW potential is concerned, A 
and Q enter uq only in the combination A Q so that a variation of the potential 
range only introduces a linear rescaling of the Q dependence of the function. 
We have checked that the graph remains qualitatively unchanged for different 
parameter settings. The first minimum of u is —0.02 for the default HCY 
potential, —0.09 for the HCY potential with eo = 0, and slightly above —0.09 for 
the SW potential. 



the result into the pde (1) yields 



df_ , rf 02 47r 2 /- 91 d K\ d 02 <f> d 2 1 
OQ 00+ Q 2 u\ Y dQe + dQuo) u dg 2 £ 



for Q away from the Q^ v Both e?oo an d d^ arc negative in the case under 
consideration [1]. 

Of the expressions appearing on the right hand side of eq. (5) the one in- 
volving the Q derivative of 1/e is of order 0(£ Seff_1 ) and so can be neglected if 
s e s < 1 as assumed. As we are looking for an effect triggered by the form of </> 
alone we do not have to consider the derivatives of the properties of the hard 
sphere reference system encoded in K, either. It is then the term involving the 
Q derivative of u that is of interest: The ideal gas contribution —l/gtoK. en- 
sures positive do2 JC so that this term is the product of duo /OQ and manifestly 
positive factors. Now assume that Qaq is less than the position of the first 
minimum of uq so that only the monotonous growth of uq towards its global 
maximum at Q = remains to be covered by the solution of the PDE: Clearly, 
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Figure 5: f(Q,g) for intermediate cutoff as a logarithmic contour plot: The data 
has been obtained for SWs with A = 3 at inverse temperature (3 = = 
0.115/e. Both the approach to the low-density boundary condition of vanishing 
/ at densities below 0.01/cr 3 and the final build-up of infinite compressibility at 
cutoffs below 10 _1 /ct have been excluded from the graph. 



as the calculation proceeds in the negative Q direction, the steeper this rise of 
uo, the more the duo/dQ term counteracts the growth of /, thereby effectively 
further delaying the onset of smoothing in Q. Most likely Qaq will thus be 
found at cutoffs so low that uq already levels off towards its limiting value of 
unity 

For the two potentials considered earlier, viz. SWs and the hcy system with 
z = 1.8/cr, fig. 4 shows that u levels off when Q (or XQ, in the case of SWs) 
is no more than about 10 _1 /(T, which is well compatible with the estimate of 
section 3.2. In addition, figs. 2 and 3 demonstrate that the transition to the 
regime where / mostly grows like l/Q, corresponding to vanishing s c s, occurs 
at similar values of the cutoff. All in all, our arguments, heuristic as they are, 
do indeed allow us to estimate Qaq in a satisfactory way. As for smoothing in 
g, on the other hand, actual numerical solution of the PDE currently is the only 
way of locating and studying Qaq{T, g). 
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5 Beyond asymptotics 



On the basis of the results presented so far one might expect numerical difficul- 
ties to first surface close to Qax, i- e., around Q ~ 10~ 1 / fJ for the potentials 
considered earlier. However, the monitoring variant of our code [9, 10, 11] that 
must be credited with first highlighting the stiffness of the equations clearly 
signals the inadequacy of the discretization grid already at much higher cutoff, 
viz., typically for 5 < Q a < 10: Indeed, the asymptotic region of large s can 
never even be reached without renouncing control of the local truncation error 
in solving the fdes, cf. section HIE of ref . [9]. 

In combination with the observed patterns of the evolution of / at inter- 
mediate and small Q illustrated in figs. 2, 3, and 5, our experience with the 
numerics of hrt leads us to propose that stiffness is not confined to that part of 
T> where the final build-up of infinite kt takes place. Indeed, fig. 5 shows that 
there are several regions of large / at higher cutoff, some of which may give rise 
to transient stiffness of the pde: Even though the analysis of part I [1] does 
not apply directly — / being bounded, asymptotic reasoning is not guaranteed 
to be valid, nor does large / imply large e any longer due to the smallness of 
Uq — , from the expressions given in part I [1] we can still deduce that do2 is 
negative and appreciable for all Q in the relevant cutoff range except very close 
to the v and that doo is likely to be rather large in modulus for / » 1 due 
to the terms linear in /. Depending on the sign of doo, large / may well prompt 
rapid further growth when Q proceeds to smaller values. Just as in section 3.3, 
such a rapid growth of / almost certainly induces an accompanying growth of 
\d 2 f/dg 2 \ on the grid, and any oscillations of the density curvature will carry 
over to df/dQ. Qualitatively the situation is then quite similar to that in the 
asymptotic region, and it seems reasonable to see this transient stiffness at in- 
termediate cutoff as preventing computations insisting on local convergence on 
a dynamically adjusted discretization mesh to ever proceed to Q <~ Qax- 

Without the backing of more formal arguments much of the above line of 
thought may seem insubstantial. There are, however, a number of numerical 
effects that provide at least indirect evidence for the point of view just laid out. 
Among those already discussed in our earlier work on hrt, the plummeting step 
sizes observed when determining the discretization grid from the local curva- 
ture of appropriate components of the solution vector [9, 11] are the most direct 
sign of stiffness at intermediate Q. Further support comes from our study of 
SWs of varying range [12]: There the peculiar shifts in the critical temperature 
whenever A is close to a simple fraction have been linked to the modulation of 
£ by the interference of c r 2 ci and (ft; and considering our remarks on the effect of 
extending the solution vector (section 3.1) it is significant that the critical point 
is accessible in a wider A range when coupling the PDE to a smaller number of 
expansion terms for taking into account the core condition, cf. section IV E of 
ref. [12] and appendix E of ref. [10]. Transient stiffness also explains why the 
lowest temperature attainable numerically, denoted 1/fcs /3max,# in refs. [10, 12], 
may well be higher than T c even though stiffness in the asymptotic region is a 
problem only for T < T c , and that the isotherms show no sign of phase separa- 
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tion for P < /3 max .# < flc, the critical temperature being known independently 
from related computations or by other methods. — There are also some more 
intricate issues related to the interplay of the ■ (where cfoo vanishes as <j) 2 ) 
with the boundaries of the cutoff ranges where the step sizes AQ are inappro- 
priate, as well as to the g dependence of the onset of smoothing in the presence 
of a local density grid refinement. Discussion of these subtle effects and their 
numerical manifestations requires a detailed presentation of appropriate meth- 
ods of data analysis on non-uniform high-resolution density grids and so falls 
outside the scope of the present report. 

6 Hard core Yukawa vs. other potentials 

In conjunction with our earlier analyses of the issues surrounding initial and 
high density boundary conditions, implementation of the core condition, and 
the peculiarities of discontinuous potentials [9, 10, 12], assertion of stiffness at 
low and intermediate Q below a certain temperature provides us with a detailed 
understanding of the numerical process of solving the hrt pde throughout V 
and has proved invaluable in interpreting numerical raw data. We close this 
short series of reports with a generally relevant sample of the kind of insight 
that can be gained on this basis, viz., a clarification of the unusually benign 
computational properties of the HCY potential: 

Throughout our numerical work we consistently found that HCY fluids of 
moderate inverse screening length like, e. g., the one with z = 1.8/u repeatedly 
used here and in ref. [9] exhibit the symptoms of stiffness only in a rather mild 
form, both for Q — > where this follows from the low value of Qaq, and at inter- 
mediate cutoff. This can be understood by noting, firstly, that the temperature 
enters the calculation only through <f> = —0w, the global normalization of which 
is accessible only in the limit Q — > 0. At any cutoff Q, the variation of <f>(k) for 
k > Q is then the only measure of the temperature available to the PDE; at the 
same time, for every one of the patches of large /, stiffness arises only below 
some characteristic temperature, coinciding with the critical one for the final 
build-up of infinite kt at Q = 0. Secondly, recalling that uo oc 4>, a look at fig. 4 
and the numbers quoted in its caption shows that the local extrema at Q > 
of (j> hcy with the default choice of e = e are substantially smaller in modulus 
than those for sws; only for much higher z do the extrema of u^ cy approach 
the SW values which they reach in the infinite- z limit. It is easily checked that 
these observations also hold in comparison with other short-ranged potentials 
like, e. g., the Lennard- Jones one: the main difference relative to sws concerns 
the phase rather than the amplitude of the oscillations. Taken together, the 
smallncss of the local extrema of <p hcy and the role the temperature plays read- 
ily explain the especially attractive numerical properties of this potential: For 
Q —* 0, the slope of mo relative to the scale set by the oscillations at higher 
Q is particularly steep, as per section 3.2 leading to especially small Qaq and 
suppression of near-discontinuities of /. At intermediate cutoff, on the other 
hand, the smallness of the local extrema vis-a-vis the global maximum at Q = 
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renders the numerics there similar to what would be seen at much higher T/T c 
in other systems, and transient stiffness poses less of a problem. At the same 
time, the z-dependence of Ug Cy immediately explains the deteriorating accuracy 
of the results for very short HCY screening length [13]. 

Support for this view comes from the numerical properties following from a 
different choice of w(r) inside the hard core, which affects the Fourier transform 
<j> and hence / and all other properties of the Q system at all cutoffs except 
in the limits Q — > oo and Q — > 0. (Independence of the final results from the 
precise choice of w is confirmed in a rather satisfactory way in some preliminary 
calculations on the Girifalco description of fullerenes [14].) The simplest such 
modification of w consists in a non-default setting of eo in eq. (2), exemplified 
by the dot-dashed curve in fig. 4 (e = 0): Just as in ref. [12], even a modest 
discontinuity of w(r) at r = a strongly affects the form of <f> and renders the 
local extrema similar to those of the sw expected, this is accompanied 

by numerical difficulties at intermediate Q similar to those discussed for SWs 
in refs. [9, 10]. In contrast, extension of the Yukawa form all the way to the 
origin — hardly unproblcmatic as it entails diverging direct correlation function 
at r — and invalidates the expansion method of taking into account the core 
condition [6] — yields the non-oscillatory form iio(Q) = z 2 /(z 2 + Q 2 ) (dotted 
curve in fig. 4) and prevents numerical solution of the fdes even at high temper- 
atures. The exceptionally attractive numerical properties of the potential (2) 
are therefore merely the result of a particular choice, shared with the original 
implementation, of w(r) inside the core and so no genuine trait of the HCY fluid. 

This finding nevertheless does not invalidate the special standing of the HCY 
system that must be taken into account in interpreting a comparison with other 
thermodynamically consistent liquid state theories on the basis of HCY results 
for 1.8/cr < z < 9/a [13]. On the other hand, the above considerations also point 
to the possibility of tuning the computational properties of some given potential 
by optimizing w(r) inside the core to reduce the local extrema at intermediate 
cutoff, an avenue largely unexplored to date the merit of which we are currently 
in no position to assess. 

The preceding clarification regarding the HCY system vis-a-vis other poten- 
tials is but one application of the detailed understanding of the hrt numerics 
presented here as well as in our earlier hrt related work. Other aspects of the 
numerics where this understanding has proved invaluable in interpreting the 
computational process and the results it yields concern the limits of the reso- 
lution in g when using extremely fine density grids, the interplay between non- 
uniform discretization grids and the location of the binodal, the local behavior 
of the solution close to the zeros i of <f>, or questions of data analysis. All in 
all, we feel that we have amassed a considerable amount of numerical experience 
and arrived at a rather detailed self-consistent perception of the computational 
process throughout all of V even below the critical temperature. Given the 
precarious nature of the HRT numerics and the not altogether unproblcmatic 
relation between the PDE and its fd approximation such an understanding is 
of prime importance if systematic mistakes are not to be introduced into the 
results unknowingly. 
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